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Abstract 


Engineering  interest  in  natural  and  ventilated  cavities  about  submerged  bodies  and  in  turbomachinery  has  led 
researchers  to  study  and  attempt  to  model  large  scale  cavitation  for  decades.  Comparatively  simple  analytical 
methods  have  been  used  widely  and  successfully  to  model  developed  cavitation,  since  the  hydrodynamics  of 
these  flows  are  often  dominated  by  irrotational  and  rotational  inviscid  effects.  However,  a range  of  more  com- 
plex physical  phenomena  are  often  associated  with  such  cavities,  including  viscous  effects,  unsteadiness, 
mass  transfer,  three-dimensionality  and  compressibility.  Though  some  of  these  complicating  physics  can  be 
accommodated  in  simpler  physical  models,  the  ongoing  maturation  and  increased  generality  of  multiphase 
Computational  Fluid  Dynamic  (CFD)  methods  has  motivated  recent  research  by  a number  of  groups  in  the 
application  of  these  methods  for  developed  cavitation  analysis.  This  paper  focuses  on  the  authors’  recent 
research  activities  in  this  area. 

The  authors  have  developed  an  implicit  algorithm  for  the  computation  of  viscous  two-phase  flows.  The  base- 
line differential  equation  system  is  the  multi-phase  Navier-Stokes  equations,  comprised  of  the  mixture  vol- 
ume, mixture  momentum  and  constituent  volume  fraction  equations.  Though  further  generalization  is 
straightforward,  a three-species  formulation  is  pursued  here,  which  separately  accounts  for  the  liquid  and 
vapor  (which  exchange  mass)  as  well  as  a non-condensable  gas  field.  The  implicit  method  developed  employs 
a dual-time,  preconditioned,  three-dimensional  algorithm,  with  multi-block  and  parallel  execution  capabili- 
ties. Time-derivative  preconditioning  is  employed  to  ensure  well-conditioned  eigenvalues,  which  is  important 
for  the  computational  efficiency  of  the  method.  Special  care  is  taken  to  ensure  that  the  resulting  eigensystem 
is  independent  of  the  density  ratio  and  the  local  volume  fraction,  which  renders  the  scheme  well-suited  to  high 
density  ratio,  largely  phase-separated  two-phase  flows  characteristic  of  developed  and  supercavitating  sys- 
tems. A dual-time  formulation  is  employed  to  accommodate  the  inherently  unsteady  physics  of  developed  and 
super-cavities.  We  have  recently  extended  the  formulation  for  compressible  constituents  to  accommodate 
analysis  of  high  speed  projectiles  and  rocket  plumes,  and  these  formulation  elements  are  also  summarized.  To 
demonstrate  the  validation  status  and  general  capabilities  of  the  scheme,  numerous  examples  are  presented. 

Nomenclature 

Symbols 


Aj 

W’  W’  W 
r*  c 

West’  ^prod 
Ci 

flux  Jacobians 
turbulence  model  constants 
mass  transfer  model  constants 
pseudo-sound  speed 

Cp 

pressure  coefficient 

cD 

drag  coefficient 

D 

source  Jacobian 

Paper  presented  at  the  RTO  A VT  Lecture  Series  on  “Supercavitating  Flows  ”,  held  at  the  von  Kdrmdn 
Institute  (VKI)  in  Brussels,  Belgium,  12-16  February  2001,  and  published  in  RTO  EN-010. 
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Subscripts,  Superscripts 

4 

i»j 

k 

L,  R 
/ 
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ng 
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+/- 
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body  diameter 

bubble  diameter 

total  energy  per  unit  volume 

flux  vectors 

frequency 

gravity  vector 

enthalpy 

source  vector 

identity  matrix 

metric  Jacobian 

transform  matrix 

turbulent  kinetic  energy 

bubble  length 

Mach  number,  similarity  transform  matrices 
mass  transfer  rates 
turbulent  kinetic  energy  production 
turbulent  Prandtl  numbers  for  k and  e 
pressure 

transport  variable  vector 

Reynolds  number 

Strouhal  number 

arc  length  along  configuration 

time  coordinate,  mean  flow  time  scale  (d/U  ) 

velocity  magnitude,  contravariant  velocity  components 

Cartesian  velocity  components 

Cartesian  coordinates 

mass  fraction 

volume  fraction,  angle-of-attack 
preconditioning  parameter 

time  derivative  preconditioning  and  transform  matrices 

turbulence  dissipation  rate,  numerical  Jacobian 

parameter,  internal  energy  per  unit  mass 

MUSCL  parameters 

eigenvalues 

molecular  viscosity 

density 

cavitation  number 
pseudo-time  coordinate 
dissipation  sensor 
curvilinear  coordinates 


single-phase  value 
coordinate  indices 

constituent  index,  pseudo-time-step  index 
dependent  variable  values  on  left  and  right  of  face 
liquid 
mixture 

non-condensable  gas 
turbulent 

condensable  vapor,  viscous 
free  stream  value 

transformed  to  curvilinear  coordinates 
production/destruction,  right/left  running 
with  respect  to  mixture 
mass  fraction  form 
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Introduction 


Multi-phase  flows  have  received  growing  research  attention  among  CFD  practitioners  due  in  large  measure  to 
the  evolving  maturity  of  single-phase  algorithms  that  have  been  adapted  to  the  increased  complexity  of  multi- 
component  systems.  However,  there  remain  a number  of  numerical  and  physical  modeling  challenges  that 
arise  in  multi-phase  CFD  analysis  beyond  those  present  in  single-phase  methods.  Principal  among  these  are 
large  constituent  density  ratios,  the  presence  of  discrete  interfaces,  significant  mass  transfer  rates,  non-equilib- 
rium interfacial  dynamics,  compressibility  effects  associated  with  the  very  low  mixture  sound  speeds  which 
can  arise,  the  presence  of  multiple  constituents  (viz.  more  than  two)  and  void  wave  propagation.  These  natu- 
rally deserve  special  attention  when  a numerical  method  is  constructed  or  adapted  for  multi-phase  flows. 

The  class  of  multiphase  flows  under  consideration  here  is  developed-  and  super-cavitating  flows,  wherein  sig- 
nificant regions  of  the  flow  are  occupied  by  gas  phase.  Depending  on  the  configuration,  such  “developed” 
[38]  cavities  are  composed  of  vapor  and/or  injected  non-condensable  gas. 

Historically,  most  efforts  to  model  large  cavities  relied  on  potential  flow  methods  applied  to  the  liquid  flow, 
while  the  bubble  shape  and  closure  conditions  were  specified.  Adaptations  of  potential  flow  methods  remain 
in  widespread  use  today  ([15]),  due  to  their  inherent  computational  efficiency,  and  their  proven  effectiveness 
in  predicting  numerous  first  order  dynamics  of  supercavitating  configurations. 

Recently,  more  general  CFD  approaches  have  been  developed  to  analyze  these  flows.  In  one  class  of  methods, 
a single  continuity  equation  is  considered  with  the  density  varying  abruptly  between  vapor  and  liquid  densi- 
ties through  an  equation  of  state.  Such  “single-continuity-equation-homogeneous”  methods  have  become 
fairly  widely  used  for  sheet  and  supercavitation  analysis  ([8],  [9],  [25],  [35],  [43],  for  example).  Although 
these  methods  can  directly  model  viscous  effects,  they  are  inherently  unable  to  distinguish  between  condens- 
able vapor  and  non-condensable  gas,  a requirement  of  ventilated  supercavitating  vehicle  analysis. 

By  solving  separate  continuity  equations  for  liquid  and  gas  phase  fields,  one  can  account  for  and  model  the 
separate  dynamics  and  thermodynamics  of  the  liquid,  condensable  vapor,  and  non-condensable  gas  fields. 
Such  multi-species  methods  are  also  termed  homogeneous  because  interfacial  dynamics  are  neglected,  that  is, 
there  is  assumed  to  be  no-slip  between  constituents  residing  in  the  same  control  volume.  A number  of 
researchers  have  adopted  this  level  of  differential  modeling,  mostly  for  the  analysis  of  natural  cavitation 
where  two  phases/constituents  are  accounted  for  ([1],  [23],  [32],  for  example).  This  is  the  level  of  modeling 
employed  here,  though  a three-species  formulation  is  used  to  account  for  two  gaseous  fields.  For  two  phases/ 
constituents  these  methods  are  very  closely  related  to  the  “single-continuity-equation-homogeneous”  methods 
addressed  above  with  interfacial  mass  transfer  modeling  supplanting  an  equation  of  state. 

Full-two-fluid  modeling,  wherein  separate  momentum  (and  in  principle  energy)  equations  are  employed  for 
the  liquid  and  vapor  constituents,  have  also  been  utilized  for  natural  cavitation  [12],  However,  in  sheet-cavity 
flows,  the  gas-liquid  interface  is  known  to  be  nearly  in  dynamic  equilibrium;  for  this  reason,  we  do  not  pursue 
a full  two-fluid  level  of  modeling. 

Sheet-  and  super-cavitating  flows  are  characterized  by  large  density  ratios  (p//pv  > 104  is  observed  in  near- 
atmospheric  water  applications),  relatively  discrete  cavity-free  stream  interfaces  and,  due  to  ventilation,  mul- 
tiple gas  phase  constituents.  Accordingly,  the  CFD  method  employed  must  accommodate  these  physics  effec- 
tively. 

Most  relevant  applications  exhibit  large  scale  unsteadiness  associated  with  re-entrant  jets,  periodic  ejection  of 
non-condensable  gas,  and  cavity  “pulsations”.  Accordingly,  we  and  others  ([23],  [26],  [32],  [35],  for  example) 
employ  a time-accurate  formulation  in  the  analysis  of  large  scale  cavitation. 
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Compressibility 

The  compressibility  of  the  liquid  and  gas  phase  constituents  can  significantly  influence  supercavitating  flow 
configurations  in  three  ways.  First,  regions  of  high  gas  volume  fraction  bubbly  flow  can  exist  if  the  cavity 
does  not  fully  envelop  the  body,  or  if  the  cavitator  is  either  poorly  designed  or  operating  off-design.  In  these 
circumstances  the  two-phase  mixture  can  have  a sound  speed  significantly  lower  than  that  of  either  liquid  or 
gas  constituents,  as  illustrated  in  Figure  1.  This  can  affect  hydrodynamic  performance  and  gives  rise  to  the 
cloud  collapse  physics  responsible  for  cavitation  damage. 


Figure  1 . Variation  of  mixture  sound  speed  in  a bubbly 
air/water  mixture.  (From  Brennan[4].) 


The  second  compressibility  issue  arises  in  candidate  rocket  propulsion  systems  for  very  high  speed  underwa- 
ter vehicles  where  exhaust  velocities  are  supersonic  (relative  to  propellant  sound  speed).  Thirdly,  very  high 
speed  underwater  projectiles  (i.e.,  unpropelled  supercavitating  darts)  can  be  deployed  at  speeds  which  are 
supersonic  with  respect  to  the  local  sound  speed  of  water.  This  application  also  requires  a compressible  CFD 
formulation  ([30],  for  example).  These  compressibility  effects  have  motivated  several  authors  to  employ  com- 
pressible formulations  ([1],  [3],  [30],  [43]  for  example)  and  the  present  authors  have  recently  extended  our 
baseline  formulation  to  accommodate  compressible  constituents. 

The  purpose  of  this  paper  is  to  present  the  numerical  methods  and  physical  models  employed,  and  many  of  the 
applications  pursued  in  the  authors’  research.  The  paper  is  organized  as  follows:  The  theoretical  formulation 
of  the  method  is  summarized,  including  the  baseline  differential  model,  inviscid  eigensystem,  physical  mod- 
els, compressibility  and  key  elements  of  the  numerical  method.  Particular  emphasis  is  placed  on  unique 
aspects  of  the  numerics  including  the  preconditioning  strategy,  resultant  eigensystem  characteristics,  and  flux 
evaluation  and  limiting  strategies  associated  with  the  resolution  of  interfaces.  This  is  followed  by  three  sets  of 
results.  The  first  set  includes  axisymmetric  steady-state  and  transient  analyses  of  natural  and  ventilated  cavita- 
tion about  several  configurations.  These  solutions  are  compared  to  experimental  measurements  to  demon- 
strate the  capability  of  the  modeling  employed.  The  second  set  of  results  includes  a variety  of  three- 
dimensional  analyses  of  sheet-  and  supercavitating  flows  of  relevance  to  high  speed  vehicles  and  turboma- 
chinery. The  third  set  of  results  includes  compressible  simulations  of  a vehicle  propulsion  plume  and  a high 
speed  supercavitating  projectile. 

Theoretical  Formulation 


Governing  Equations 

The  baseline  governing  differential  system  employed  is  cast  in  Cartesian  coordinates  as: 
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where  oq  and  ant!  represent  the  liquid  phase  and  non-condensable  gas  volume  fractions,  and  mixture  density 
and  mixture  turbulent  viscosity  are  defined  as: 


Pm  = P/a/  + Pvav  + Pnga„g 

PmCgk2 


M- 


m,  t 


(2) 


In  this  baseline  formulation,  the  density  of  each  constituent  is  taken  as  constant.  The  mass  transfer  rates  from 
vapor  to  liquid  and  from  liquid  to  vapor  are  denoted  m and  m , respectively.  Mass  transfer  terms  appear  in 
the  mixture  continuity  equation  because  this  equation  is  a statement  of  mixture  volume  conservation.  Also, 
note  that  each  of  the  equations  contains  two  sets  of  time-derivatives  - those  written  in  terms  of  the  variable  “t” 
correspond  to  physical  time  terms,  while  those  written  in  terms  of  “x”  correspond  to  pseudo-time  terms  that 
are  employed  in  the  time-iterative  solution  procedure.  The  forms  of  the  pseudo-time  terms  will  be  discussed 
presently. 


In  the  development  of  the  differential  system  presented  above,  a number  of  physical,  numerical,  and  practical 
issues  were  considered.  First,  a mixture  volume  continuity  equation  is  employed  rather  than  a mixture  mass 
equation.  This  initial  choice  was  made  based  on  the  authors’  experience  that  the  nonlinear  performance  of 
segregated  pressure  based  algorithms  [17]  is  improved  by  doing  so  for  high  density  ratio  multi-phase  systems. 
Because  of  this  choice,  neither  a physical  time  derivative  nor  mixture  density  appears  in  the  continuity  equa- 
tion, although  the  mixture  density  can  vary  in  space  and  time.  To  render  the  system  hyperbolic  and  to  facilitate 
the  use  of  time-marching  procedures,  we  then  introduce  a pseudo-time  derivative  term  (signified  by  “x”)  in  the 
mixture  continuity  equation,  a strategy  that  derives  from  the  work  of  Chorin  [7]  and  others. 


Second,  corresponding  artificial  time-derivative  terms  are  also  introduced  in  the  component  phasic  continuity 
equations,  which  ensures  that  the  proper  differential  equation  (in  non-conservative  form)  is  satisfied.  That  is, 
combining  equations  [l]a  and  [l]c: 
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Inclusion  of  such  ‘'phasic  continuity  enforcing”  terms  has  a favorable  impact  on  the  nonlinear  performance  of 
multi-phase  algorithms  when  mass  transfer  is  present  [34], 

Third,  we  desired  an  eigensystem  that  is  independent  of  density  ratio  and  volume  fractions  so  that  the  perfor- 
mance of  the  algorithm  would  be  commensurate  with  that  of  single-phase  for  a wide  range  of  multi-phase 
conditions.  These  considerations  give  rise  to  the  preconditioned  system  in  equation  [1], 

In  generalized  coordinates,  equations  [1]  can  be  written  in  vector  form  as: 


0, 


where  the  primitive  solution  variable,  flux,  and  source  vectors  are  written: 
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and  J is  the  metric  Jacobian,  J = d(x,  y,  z)/d(£,  rj,  Q . 
Matrix  Te  is  defined  by: 
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and  the  preconditioning  matrix,  T,  takes  the  form: 
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where  APl  = P/-pv  and  Ap2  = png-pv.  The  compatibility  condition,  oq  + av  + ang  = 1 , is  incorporated 
implicitly  in  definitions  6 and  7. 
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Eigensystem 

Of  interest  in  the  construction  and  analysis  of  a scheme  to  discretize  and  solve  equation  [4]  is  its  inviscid 
eigensystem.  In  particular,  the  eigenvalues  and  eigenvectors  of  matrix  Aj  are  required,  where 


Aj  = r'Aj 


A:=^. 

3Q 


(8) 


Aj,  r"1 , and  Aj  can  be  computed  straightforwardly,  and  expressions  for  these  are  available  in  [18].  The 
eigenvalues  and  eigenvectors  of  Aj  can  be  found  by  first  considering  the  reduction  of  equation  [4]  to  a single- 
phase system.  With  a/  = 1,  pm  = p/  = pv  = png  (=  constant),  m =0,  equation  [4]  collapses  to  the  widely  used 
single-phase  “pseudo-compressibility”  scheme,  which  can  be  written  for  inviscid  flow  as 
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Flux  Jacobian  matrix  Aj  0 has  a well  known  form  ([28],  for  example)  and  is  also  given  in  [18]. 
Comparing  the  expressions  for  A^  and  Aj , one  can  write: 
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Diagonalizing  A- : 


Aj  = MjAjMj”1  . (11) 

The  elements  of  the  diagonal  matrix  Aj  are  the  eigenvalues  of  Aj . Similarity  transform  matrices  Mj  and  Mj  1 
contain  the  right  and  left  eigenvectors  of  Aj . Forms  for  these  matrices  are  sought.  Using  equation  [10],  we  can 
write  equation  [ 1 1 ] as 
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(12) 


This  brief  analysis  illustrates  that  the  inviscid  eigenvalues  of  the  present  preconditioned  multi-phase  system 
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are  the  same  as  the  standard  single-phase  pseudo-compressibility  system  with  two  additional  eigenvalues 
introduced,  IT,  Uj.  That  is. 


Ai-(Ui.Ui.Uj+Cj.Uj-Cj.Uj.U))T 
C,  ■ >?+s2«iAi) 

Equation  [12]  also  illustrates  that  a complete  set  of  linearly  independent  eigenvectors  exists  for  the  present 
three-component  system.  These  results  generalize  for  an  arbitrary  number  of  constituents.  The  eigenvalues  are 
seen  to  be  independent  of  the  volume  fractions  and  density  ratio.  This  is  not  the  case  for  other  choices  of  pre- 
conditioning matrix,  T,  or  if  a mixture  mass  conservation  equation  is  chosen  instead  of  the  mixture  volume 
equation. 

The  local  time-steps  and  matrix  dissipation  operators  presented  below  are  derived  from  the  inviscid  multi- 
phase eigensystem  above,  which  has  been  shown  to  be  closely  related  to  the  known  single-phase  eigensystem. 
This  has  had  the  practical  advantage  of  making  the  single-phase  predecessor  code  easier  to  adapt  to  the  multi- 
phase system. 

Physical  Modeling 
Mass  Transfer 

For  transformation  of  liquid  to  vapor,  m is  modeled  as  being  proportional  to  the  liquid  volume  fraction  and 
the  amount  by  which  the  pressure  is  below  the  vapor  pressure.  This  model  is  similar  to  that  used  by  Merkle  et 
al.  [23]  for  both  evaporation  and  condensation.  For  transformation  of  vapor  to  liquid,  m , a simplified  form  of 
the  Ginzburg-Landau  potential  is  employed: 


• - cdestPva/MIN[0,p-pv] 

m - - 

(1/2P/U  l)toe 

m+  = CprodPv(a/— angH]  -a/~Kng) 


(14) 


In  this  work,  Cdest  and  Cpt.0(j  are  empirical  constants  (here  Ct|est  = 1 05,  Cprod  = 1 05).  an„  appears  in  the  produc- 
tion term  to  enforce  that  m — » 0 as  av  — > 0 . Both  mass  transfer  rates  are  non-dimensionalized  with  respect  to 
a mean  flow  time  scale. 


The  mass  transfer  model  presented  in  equation  [14]  retains  the  physically  observed  characteristic  that  cavity 
sizes,  and  thereby  the  dynamics  of  the  two- fluid  motion,  are  nearly  independent  of  liquid-vapor  density  ratio. 
One  outcome  of  the  eigensystem  characteristics  summarized  above  is  that  the  numerical  behavior  of  the  code 
parallels  this  physical  behavior.  This  is  illustrated  in  Figure  2 There,  the  convergence  histories  are  provided 
for  six  simulations  of  cavitating  flow  over  a hemispherical  forebody  with  cylindrical  afterbody  (1/2  caliber 
ogive,  o = 0.3,  from  series  of  results  presented  below).  As  the  density  ratio  is  increased  from  1 to  10,  the  con- 
vergence history  is  modified  somewhat,  but  beyond  p/pv  = 10,  the  performance  of  the  solver  is  virtually  inde- 
pendent of  density  ratio  up  to  p//pv  = 105.  This  behavior  is  consistent  with  the  modeled  physics  of  this 
problem  as  shown  in  Figure  3.  There,  drag  coefficient  and  number  of  predicted  vaporous  cells  are  seen  to 
reach  a nearly  constant  value  at  p//pv  > 10. 

As  the  authors  have  evolved  this  work  into  the  analysis  of  turbomachinery,  it  has  become  evident  that  the 
present  mass  transfer  model  has  some  shortcomings  in  that  it  does  not  accommodate  thermal  effects  on  cavita- 
tion breakdown  performance  in  turbomachinery  [5],  This  is  discussed  further  in  the  pump  cavitation  results 
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section  below. 


Figure  2.  Comparison  of  convergence  histories  with  den-  FlSure  1 Comparison  of  predicted  drag  coefficient 

sity  ratio  for  hemispherical  forebody  simulation.  and  number  of  vaporous  cells  with  density  ratio  for 

hemispherical  forebody  simulation. 


Turbulence  Closure 

A high  Reynolds  number  form  k-e  model  (Jones  and  Launder  [16])  with  standard  wall  functions  is  imple- 
mented to  provide  turbulence  closure: 


s<p”k)+4;<p"'kni) 


9 f bm.  idk  3 
9xjVPrtk9xjJ 


+ P-  pe 


(15) 


3 ( 

3T(P' 


'e)  + S-(p 


9 fbm.  t9e 


m£lV  9x1  Pr 


te0'M 


y + [C|P-C2pe,(5) 


As  with  velocity,  the  turbulence  scalars  are  interpreted  as  being  mixture  quantities.  Other  two-equation  turbu- 
lence models  have  also  been  employed  as  discussed  below. 


Compressible  Constituent  Phases 


As  mentioned  above,  a complicating  phenomena  associated  with  underwater  multiphase  flows  is  the  presence 
and  effect  of  compressibility  in  a flow  that  is  largely  incompressible.  To  directly  model  flows  containing 
homogeneously  mixed  bubbly  flow  regions  and/or  shock-expansion  wave  systems,  a suitable  representation 
of  compressible  flow  is  necessary.  Also,  in  flows  of  relevance  here,  liquid  vapor  mass  transfer  is  important 
and  by  virtue  of  ventilation,  a noncondensable  gas  phase  is  present. 


For  the  purposes  of  analysis  and  development  of  appropriate  preconditioning,  the  inviscid  one-dimensional 
two-phase  form  of  the  governing  equations  are  presented  in  equation  [16].  Here,  with  no  loss  of  generality,  the 
equations  are  presented  in  mass-fraction  form.  In  this  form,  the  governing  equations  for  two-phase  flow 
resemble  the  equations  employed  in  single-phase  multi-component  reacting-gas-mixture  flows.  The  eigenval- 
ues of  the  associated  system  and  the  preconditioning  fonns  are  well  known  and  have  been  widely  used  [40], 
The  equations  corresponding  to  a preferred  set  of  primitive  variables  may  be  achieved  and  then  interpreted  in 
the  context  of  multi-phase  flows  of  interest.  A more  complete  development  with  applications  is  given  in  refer- 
ences [20],  [21]  and  [41], 
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Compressible  Equations  of  Motion 

The  fully  compressible  form  of  the  two-phase  equations  in  ID  is  given  as, 


dQy  , dE 
dr  dx 


- 0 


which  represents  the  following  set  of  equations  in  conservation  form, 


(16) 


dPYv  dpYvu  n 
dx  dx 

(17) 

dp  Y,  dp  Y,u 
dr  + dx 

(18) 

2 

dp  u dp  if  dp  _ 
dr  dx  dx 

(19) 

de  , d(e  + p)u  _ A 
dr  dx 

(20) 

The  mixture  density  is  defined  as: 

P = Pv  + P/  = p,a/  + pvav  (21) 

where, 

Pv  = pYv  = Pvav  and  P/  = py,  = p,a,  (22) 

The  total  internal  energy  of  the  two-phase  mixture  may  be  expressed  as, 

e = pe  + ipu2  = ph-p  + ipu2  (23) 

and  the  specific  enthalpy  of  the  mixture  is  given  as, 

h=  I hjYj  (24) 

i = /,  v 

and  h[  and  hy  are  the  phasic  enthalpies,  which  are  also  generally  known  functions  of  the  temperature  and  the 
pressure.  The  p indicates  where,  in  equations  [21]  and  [22],  the  mass  fraction  based,  Dalton  model,  of  species 
density  is  to  be  used,  rather  than  the  volume  fraction  based,  Amagat  model.  For  development  of  the  compress- 
ible preconditioning  form,  the  Dalton  model  for  intensive  properties  is  used  here. 

The  system  in  equations  [17]-[20]  may  equivalently  be  expressed  in  the  following  vector  form: 

fe^  + 9E  = 0 (25) 

or  dx 

where  we  have  adopted  a set  of  primitive  variables,  Q , as  the  primary  dependent  vector: 
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The  corresponding  system  flux  Jacobian  is  defined  as: 


A = 


uY  — 
vDp 

PU  + UYv^“ 

PYV 

uY  — 
VDT 

Q 

1 

Q 

6 

uV^P 

uY/Dp 

-Pu  + uY/^- 

PY/ 

uY 

Y/dt 

Q 

Q 

Q 
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DT 

Q 

Q 

Q 

pu 


Dh 

Dp 
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Dh 

DY. 
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°DY, 


,.  , Dh 

p(h0  + u2)  pu^ 


+ uh 


Dp 

°DT 


(26) 


Note  that  a standard  notation  has  been  used  to  indicate  partial  differentiation  with  respect  to  the  variables  in 
Q . The  eigenvalues  of  the  above  system  are: 


A.([Fe]  A)  = u,  u,  u+c 


(27) 


where  the  speed  of  sound  is  given  by:. 


,2  = 


A 

Dp 


Dh 

PDT 

Q 

Dh 

fi  dh 

l-o  — 

\ 

.DT 

_ DT 

PDp 

Q 

Q 

Q 

V 

Q 

For  the  multi-phase  system,  the  above  partial  derivatives  need  to  be  obtained  in  terms  of  the  known  volume 
fraction  based  properties,  i.e.,  py  = p (p,  T)  and  p,  = p/(p,  T) . Their  evaluation  is  therefore  not  as  straight- 
forward as  for  single-phase  multi-component  mixtures. 

The  relation  for  the  isothermal  sound  speed  is  given  in  equation  [28]. 
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Preconditioned  Compressible  Equations  of  Motion 


The  preconditioned  version  of  equation  [25]  may  be  written  as: 


fdQ  OE 
dr  dx 


= 0 


where  the  preconditioning  matrix  is  given  by: 


(28) 


(29) 
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The  eigenvalues  of  the  preconditioned  system  can  be  shown  to  be: 


(30) 


X(F'a)  = u,u,i 


<'+(l)>4-(92)2+- 


where  the  preconditioned  pseudo-sound-speed  is: 


(o')2  = 


dp' 


dh 

dT 


dh 

PdT 


dT 


1 -p 


dh 

dp 


Definition  of  Pseudo-Sound  Speed 


(31) 


(32) 


The  definition  of  the  pseudo-sound  speed  is  selected  to  render  the  eigenvalues  of  the  preconditioned  system 
well  conditioned.  Examining  equation  [31],  we  note  that  this  is  readily  achieved  by  choosing: 


(c')2  = Min(V2,  c2)  (33) 

where  V is  local  convective  velocity  magnitude.  This  form  corresponds  to  the  standard  inviscid  choice  of  this 
parameter.  From  equation  [32],  we  may  thus  express  the  preconditioning  parameter  in  equation  [30]: 
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Numerical  Method 

The  baseline  numerical  method  is  evolved  from  the  UNCLE  code  of  Taylor  and  his  co-workers  at  Mississippi 
State  University  ([39],  for  example).  UNCLE  is  based  on  a single-phase,  pseudo-compressibility  formulation. 
Roe-based  flux  difference  splitting  is  utilized  for  convection  term  discretization.  An  implicit  procedure  is 
adopted  with  inviscid  and  viscous  flux  Jacobians  approximated  numerically.  A block-symmetric  Gauss-Seidel 
iteration  is  used  to  solve  the  approximate  Newton  system  at  each  time-step. 


The  multi-phase  extension  of  the  code  retains  these  underlying  numerics  but  also  incorporates  two  volume 
fraction  transport  equations,  mass  transfer,  non-diagonal  preconditioning,  flux  limiting,  dual-time-stepping, 
and  two-equation  turbulence  modeling.  The  fully  compressible  version  incorporates  an  energy  conservation 
equation. 


Discretization 


The  transformed  system  of  governing  equations  is  discretized  using  a cell  centered  finite  volume  procedure. 
Flux  derivatives  are  computed  as 


dE/d£,  = (Ej+i/2-  Ei-1/2)  , (35) 

with  similar  expressions  for  dF/dr) , dG/d £ , and  the  corresponding  viscous  fluxes.  The  inviscid  numerical 
fluxes  are  evaluated  using  a flux  difference  splitting  procedure  [45]: 


Ei+i/2=  J[E(Q!;1/2)  + E(Qj.I/2)  + 


A! 


I(Qh 


(36) 


i+l/2’  Qi+1/2)  * (Qi+1/2  <W] 


where,  with  the  non-diagonal  preconditioner  used  here,  the  matrix  dissipation  operator  is  defined  by 
|a|=  f(M|A|M“'). 


R L 

The  extrapolated  Riemann  variables,  Qi+1/2  and  Qj+|/2  are  obtained  using  a MUSCL  procedure  ([2],  for 
example): 


Qw/2  = Qi+1  - 1[(  l-K)(Qi+2  - Qi+1 ) + ( l+K)(Qi+i  - Qi)] 

Qi+1/2  = Qi  +|[(l-K)(Qi  -Qj.,  ) + (l+K)(Qi+]  -Qj)] 

(37) 

For  first  order  accuracy,  4>  = 0.  The  choice  d = 1,  k = 1/3,  yields  the  third  order  accurate  upwind  bias  scheme 
used  for  the  results  presented  in  this  paper. 

The  flows  of  interest  here  typically  contain  regions  with  sharp  interfaces  between  liquid  and  gas  phases.  In 
addition,  compressible  flows  admit  shock  waves  and  contact  discontinuities.  Accordingly,  higher  order  dis- 
cretization practices  are  required  to  retain  adequate  interface  fidelity  in  the  simulations.  This  is  particularly 
important  in  three-dimensional  super-cavitating  vehicle  or  control  surface  computations  such  as  those  pre- 
sented below.  There,  predicted  lift  and  drag  can  be  severely  over-predicted  if  liquid  phase  (and  its  much 
higher  inertia)  diffuses  numerically  into  low-lift  gaseous  regions  of  the  lifting  surface. 
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Attendant  to  the  third  order  upwind  bias  scheme  employed  are  overshoots  in  solution  variables  at  these  inter- 
faces. These  can  be  highly  destabilizing,  particularly  for  the  volume  fraction  equations,  if  sufficient  mass 
transfer  or  non-condensable  vapor  is  present  to  yield  a;  -» 0 locally.  To  ameliorate  this  difficulty,  the  flux 
evaluation  is  rendered  locally  first  order  in  the  presence  of  large  gradients  in  a;  or  o^g.  In  solutions  obtained 
for  incompressible  phasic  constituents,  this  is  affected  through  the  use  of  a “dissipation  sensor”  in  the  spirit  of 
Jameson  et  al.  [14],  Specifically,  a sensor  is  formulated  for  each  coordinate  direction  as 

aj+1  - 2aj  + oCj., 
aj+1  + 2aj  + a-  , 

This  parameter  is  very  small  except  in  the  immediate  vicinity  of  liquid-gas  interfaces.  In  this  work,  the  higher 
order  component  of  the  numerical  flux  in  equation  [37],  i.e.  term  <]),  is  multiplied  by  (1-Vj).  In  compressible 
solutions,  shocks  and  contact  discontinuities  are  admissible.  Hence,  it  is  necessary  to  limit  interpolations 
based  on  the  entire  primitive  vector.  A standard  form  of  the  van  Albada  limiter  [13]  has  been  employed. 

To  illustrate  the  foregoing  discretization  issues,  consider  adjacent  parallel  streams  of  two  constituents.  In  the 
absence  of  shear,  if  a flow-aligned  Cartesian  mesh  is  employed,  the  mixing  layer  interface  will  be  perfectly 
preserved  using  the  present  modeling  (no  mass  diffusion).  This  is  true  independent  of  the  density  ratio  of  the 
streams  and  whether  first  or  higher  order  discretization  is  employed.  However,  if  the  interface  encounters  a 
region  of  significant  grid  non-orthogonality,  the  interface  will  be  smeared. 

We  consider  a two-dimensional  grid  slice  extracted  from  the  three-dimensional  fin  computation  presented 
below,  as  illustrated  in  Figure  4a.  Two-dimensional  inviscid  computations  of  a high  density  ratio  mixing 
stream  were  computed  on  this  grid.  Gas  and  liquid  (p/p„  = 1000)  were  injected  axially,  at  the  same  velocity, 
above  and  below  an  inlet  location  seen  in  Figure  4b.  No  solid  boundaries  are  specified,  so  ideally,  the  interface 
would  be  perfectly  preserved  through  the  domain.  However,  the  initially  sharp  mixing  layer  interface  encoun- 
ters severe  grid  non-orthogonality,  as  it  does  in  full  scale  three-dimensional  computations,  even  well  away 
from  solid  boundaries,  due  to  grid  topology.  Though  this  case  is  particularly  pathological,  regions  where  local 
grid  quality  suffers  due  to  the  geometry  of  the  problem  are  inescapable  in  structured  multiblock  analyses  of 
complex  aerodynamic  configurations. 

Downstream  of  the  grid  “stripe”  associated  with  the  fin  leading  edge,  Figures  4b  and  c show  that  the  interface 
is  badly  smeared  when  first  order  discretization  is  employed  (0.05  < a/<  0.95  for  11-12  nodes).  This  is  largely 
remedied  when  the  3rd  order  flux  difference  with  dissipation  sensor  is  employed  (0.05  < a / < 0.95  for  3-4 
nodes). 

The  dissipation  sensor  used  does  not  completely  eliminate  overshoots,  as  illustrated  in  Figure  4c.  If  the  vol- 
ume fractions  are  “clipped”  at  1.0  after  each  pseudo-time-step  the  impact  on  the  solution  accuracy  is  minimal, 
as  illustrated  in  Figure  4c.  But,  clipping  causes  the  non-linear  convergence  to  flat-line,  so  in  general,  we 
accept  the  modest  overshoots  in  volume  fraction. 

Second  order  accurate  backward  differencing  is  used  to  discretize  the  physical  transient  term  as 

3Q  (3Qn+l'k- 4Qn  + Qn  ') 
e at  e 2At 

where  index  n designates  the  physical  time-step.  Second  order  accurate  central  differencing  is  utilized  for  the 
viscous  flux  terms. 

Implicit  Solution  Procedure 

Adopting  Euler  implicit  differencing  for  the  pseudo-transient  term,  equation  [4]  can  be  written  in  A-form  as 


(39) 


(38) 
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(40) 


where  AQ  = Q - Q , with  index  n designating  the  physical  time-step  and  index  k designating  the 
pseudo-time-step.  Terms  Aj , Aj  and  D are  the  inviscid  flux,  viscous  flux,  and  source  Jacobians.  The  source 
Jacobian  is  evaluated  analytically  as  illustrated  below.  The  inviscid  and  viscous  flux  Jacobians  are  evaluated 
numerically  as: 


(A;AQ)j+]/2  = AQl  + -t— -AQr  « 

3Ql  OQr 


dE 


AQ 


+ -^P-AQi+i 


3Ql  3Qr  (41) 

~E(Ql  + e,  Qr)  - E(Ql,  Qr)JAq.  + 

~E(Ql,  Qr  + e)  - E(Ql,  Qr) jAQ. : ) 

where  e is  taken  as  the  square  root  of  a floating  point  number  close  to  the  smallest  resolvable  by  the  hardware. 


Figure  4.  a)  Two-dimensional  grid  slice  extracted  from 
three-dimensional  control  fin  model,  b)  Predicted  liquid 
volume  fraction  contours  using  present  third  order  dis- 
cretization. c)  Predicted  liquid  volume  fraction  profiles 
for  various  grids  and  discretizations. 
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The  pseudo-time-step  is  defined  based  on  the  spectral  radii  of  U1  Aj  as 

Ax  = CFL  . (42) 

xb+c,i 

j 

For  steady  state  computations,  the  physical  time-step,  At,  is  set  to  infinity,  and  a CFL  number  of  3-5  is  typi- 
cally used. 

In  the  present  work,  we  seek  to  resolve  transient  features  characterized  by  the  quasi-periodic  shedding  of  vor- 
ticity  and  gas  from  the  cavity.  Strouhal  numbers  on  the  order  of  0.1  are  encountered  in  practice.  Accordingly, 
non-dimensional  physical  time-steps  on  the  order  of  .005  - .01  are  required.  A CFL  number  of  3-5  is  typically 
used  for  the  inner  iterates  in  transient  computations.  As  illustrated  below,  this  gives  rise  to  a dual-time  scheme 
that  provides  a 1 to  3 order-of-magnitude  drop  in  residuals  in  5-10  pseudo-time-steps. 

Optimum  non-linear  convergence  is  obtained  using  a pseudo-compressibility  parameter,  p /U  M = 10. 

Upon  application  of  the  discretization  and  numerical  linearization  strategies  defined,  equation  [40]  represents 
an  algebraic  system  of  equations  for  AQ.  This  block  (6x6  blocks)  septadiagonal  system  is  solved  iteratively 
using  a block  symmetric  Gauss-Seidel  method.  Five  sweeps  of  the  BSGS  scheme  are  applied  at  each  pseudo- 
time-step. 

Source  Terms 


Following  the  strategy  of  Venkateswaran  et  al.  [40]  for  the  numerical  treatment  of  mass  transfer  source  terms 
in  reacting  flow  computations,  we  identify  a source  and  sink  component  of  the  mass  transfer  model  and  treat 
the  sink  term  implicitly  and  the  source  term  explicitly.  Specifically,  with  reference  to  equation  [1]  we  have 


H 


+/- 


1 \ . +/- 
— , 0,  0,  0,  m 

pj 


With  m from  equation  [14]  we  have 


(43) 


AH"  = D'AQ,  with  D"  = OH'/OQ  . (44) 

It  can  easily  be  shown  that  the  non- zero  eigenvalue  of  D"  is 

VD-)  = Cdt„p>[(I-i)a,  + (I),p-pl)]  (45) 

which  is  less  than  zero  for  pv  < pj.  Hence,  the  identification  of  m as  a sink  is  numerically  valid.  Implicit  treat- 
ment of  this  term  provides  that  a/  approaches  zero  exponentially,  so  that  cases  like  those  considered  below, 
where  significant  mass  transfer  results  in  extremely  low  liquid  volume  fractions,  remain  stable. 

The  production  mass  transfer  term  is  treated  explicitly.  A relaxation  factor  of  0.1  is  applied  at  each  pseudo- 
time-step  to  keep  this  term  from  destabilizing  the  code  in  early  iterations. 

Turbulence  Model  Implementation 

The  turbulence  transport  equations  are  solved  subsequent  to  the  mean  flow  equations  at  each  pseudo-time- 
step.  A first  order  accurate  flux  difference  splitting  procedure  similar  to  that  outlined  above  for  the  mean  flow 
equations  is  utilized  for  convection  term  discretization.  The  k and  e equations  are  solved  implicitly  using  con- 
ventional implicit  source  term  treatments  and  a 2x2  block  symmetric  Gauss-Seidel  procedure. 
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Boundary  Conditions 

Velocity  components,  volume  fractions,  turbulence  intensity,  and  turbulence  length  scale  are  specified  at 
inflow  boundaries  and  extrapolated  at  outflow  boundaries.  Pressure  distribution  is  specified  at  outflow  bound- 
aries (p=0  for  single-phase  or  non-buoyant  multi-phase  computations)  and  extrapolated  at  inflow  boundaries. 
At  walls,  pressure  and  volume  fractions  are  extrapolated,  and  velocity  components  and  turbulence  quantities 
are  enforced  using  conventional  wall  functions.  Boundary  conditions  are  imposed  in  a purely  explicit  fashion 
by  loading  “dummy”  cells  with  appropriate  values  at  each  pseudo-time-step. 

Parallel  Implementation 

The  multiblock  code  is  instrumented  with  MPI  for  parallel  execution  based  on  domain  decomposition.  Inter- 
block communication  is  affected  at  the  non-linear  level  through  boundary  condition  updates  and  at  the  linear 
solver  level  by  loading  AQ  from  adjacent  blocks  into  “dummy”  cells  at  each  SGS  sweep.  This  is  not  as 
implicit  as  solving  an  optimally  ordered  linear  system  for  the  entire  domain  at  each  SGS  sweep,  but,  as  dem- 
onstrated below,  this  potential  shortcoming  will  not  deteriorate  the  non-linear  performance  of  the  scheme  if 
the  linear  solver  residuals  are  reduced  adequately  at  each  pseudo-time-step.  The  authors  routinely  employ  48- 
80  processors  for  three-dimensional,  unsteady  analyses. 

Results 

Three  sets  of  results  from  the  authors’  recent  research  are  presented.  The  first  set  includes  axisymmetric 
steady-state  and  transient  analyses  of  natural  and  ventilated  cavitation  about  several  configurations.  These 
solutions  are  compared  to  experimental  measurements  to  demonstrate  the  capability  of  the  modeling 
employed.  The  second  set  of  results  includes  a variety  three-dimensional  analyses  of  sheet-  and  supercavitat- 
ing  flows  of  relevance  to  high  speed  vehicles  and  turbomachinery.  The  third  set  of  results  includes  compress- 
ible simulations  of  a vehicle  propulsion  plume  and  a high  speed  supercavitating  projectile. 

Axisymmetric  Analyses 

Steady  State  and  Transient  Natural  and  Ventilated  Cavitation  about  a Series  of  Axisymmetric 
Forebodies 

Due  to  the  authors’  principal  research  interests  in  supercavitating  vehicles,  an  emphasis  of  our  CFD  develop- 
ment and  validation  efforts  has  been  on  axisymmetric  bodies.  In  this  regard,  we  have  used  data  due  to  Rouse 
and  McNown  [29]  as  well  as  that  due  to  Stinebring  et  al.  [36],  [37]  for  validation  purposes. 

Rouse  and  McNown  [29]  carried  out  a series  of  experiments  wherein  cavitation  induced  by  convex  curvature 
aft  of  various  axisymmetric  forebodies  with  cylindrical  afterbodies  was  investigated.  At  low  cavitation  num- 
bers, these  flows  exhibit  natural  cavitation  initiating  near  or  just  aft  of  the  intersection  between  the  forebody, 
or  cavitator,  and  the  cylindrical  body.  For  each  configuration,  measurements  were  made  across  a range  of  cav- 
itation numbers,  including  a single  phase  case  (large  a).  Surface  static  pressure  measurements  were  taken 
along  the  cavitator  and  after-body.  Photographs  were  also  taken  from  which  approximate  bubble  size  and 
shape  were  deduced. 

Several  of  the  Rouse-McNown  configurations  were  analyzed.  These  included  0-caliber  (blunt),  1/4-caliber, 
1/2-caliber  (hemispherical),  1-caliber,  and  2-caliber  ogives  and  conical  (22.5°  cone  half-angle)  cavitator 
shapes.  The  experiments  were  performed  at  Reynolds  numbers  greater  than  100,000  based  on  maximum  cavi- 
tator (i.e.,  after-body)  diameter.  A value  of  Re  = 136000  was  used  for  the  simulations.  In  order  to  properly 
assess  grid  resolution  requirements,  a range  of  grid  sizes  was  used.  For  the  hemispherical  and  conical  config- 
urations, grid  sizes  of  65x17,  129x33  and  257x65  were  run.  Figure  5 demonstrates  that  differences  between 
predicted  surface  pressures  for  the  medium  and  fine  meshes  are  small.  The  fine  meshes  were  used  for  all  sub- 
sequent calculations  presented  here.  For  the  blunt  fore-body,  a two-block  grid  topology  had  to  be  used,  and  a 
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mesh  consistent  with  the  resolution  and  clustering  of  the  other  head-forms  was  utilized  (65x49,  257x65  for 
blocks  1 and  2;  see  Figure  9b). 


Figure  5.  Comparison  of  predicted  surface  pressure  distributions  for 
naturally  cavitating  axisymmetric  flow  over  a conical  cavitator-cylindrical 
after-body  configuration,  c = 0.3.  Coarse  (65x17),  medium  (129x33)  and 
fine  (257x65)  mesh  solutions  are  plotted. 

Figures  6 through  24  show  sample  results  for  these  axisymmetric  computations.  Figure  6a  shows  predicted 
and  measured  surface  pressure  distributions  at  several  cavitation  numbers  for  the  1 -caliber  ogive  forebody 
with  cylindrical  afterbody.  As  the  cavitation  number  is  decreased  from  near  a critical  “inception”  value,  a cav- 
itation bubble  forms  and  grows.  The  presence  of  the  bubble  manifests  itself  as  a decrease  in  magnitude,  flat- 
tening and  lengthening  of  the  pressure  minimum  along  the  surface.  Also,  bubble  closure  gives  rise  to  an 
overshoot  in  pressure  recovery  due  to  the  local  stagnation  associated  with  free-stream  liquid  flowing  over  the 
convex  curvature  at  the  aft  end  of  the  bubble.  The  code  is  seen  to  accurately  capture  these  physics  as  evi- 
denced by  the  close  correspondence  between  predicted  and  measured  pressure  distributions.  Figure  6b  illus- 
trates the  qualitative  physics  as  captured  by  the  model.  There,  surface  pressure  contours,  field  liquid  volume 
fraction  contours,  selected  streamlines  and  the  grid  used  are  shown  for  the  1 -caliber  case  at  c = 0. 15.  For  this 
case,  the  cavitation  bubble  is  quite  long  (L/d  > 3).  As  with  all  large  cavitation  bubbles,  the  closure  region  is 
characterized  by  an  unsteady  “re-entrant”  jet.  The  significant  flow  recirculation  and  associated  shedding  of 
vorticity  and  vapor  in  these  flows  require  that  a transient  simulation  be  carried  out.  This  is  discussed  further 
below.  The  solution  depicted  in  Figure  6b  represents  a snapshot  in  time  of  an  unsteady  simulation. 

Figures  7-9  provide  similar  comparisons  for  hemispherical,  conical  and  blunt  ogives.  Specifically,  Figures  7a, 
8a  and  9a  show  comparisons  between  predicted  and  measured  surface  pressure  distributions  for  these  three 
configurations  at  a range  of  cavitation  numbers.  In  Figures  7b,  8b  and  9b,  surface  pressure  contours,  field  liq- 
uid volume  fraction  contours,  selected  streamlines  and  the  grid  are  illustrated  for  selected  cavitation  numbers. 
In  each  case,  a discrete  bubble  shape  is  observed,  but  the  aft  end  of  the  predicted  bubble  does  not  exhibit  a 
smooth  “ellipsoidal”  closure.  Indeed,  due  to  local  flow  reversal  (reentrant  jet),  liquid  is  swept  back  underneath 
the  vapor  pocket.  The  pressure  in  this  region  retains  the  nearly  constant  free  stream  liquid  flow  value 
impressed  through  the  bubble.  Similar  closure  region  observations  have  been  made  by  Chen  and  Heister  [6] 
and  others. 
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Figure  6a.  Comparison  of  predicted  and  measured  sur-  Figure  6b.  Predicted  liquid  volume  fraction  and 

face  pressure  distributions  at  several  cavitation  num-  surface  pressure  contours,  selected  streamlines  and 

bers  for  a 1 -caliber  ogive  forebody.  computational  grid  for  a 1 -caliber  ogive  forebody, 


As  with  the  1 -caliber  ogive,  the  predicted  pressure  distribution  for  the  hemispherical  cavitator  is  very  well 
predicted,  but  Figure  8b  illustrates  that  the  simulations  underpredict  the  length  of  the  bubble  for  the  22.5° 
half-angle  conical  head  form  at  all  cavitation  numbers,  though  qualitative  trends  remain  correctly  predicted. 
Although  fairly  good  quantitative  agreement  is  obtained  for  the  blunt  fore-body  at  low  cavitation  numbers, 
shortcomings  of  the  present  modeling  are  evident.  First,  even  the  single  phase  pressure  distribution  shows  sig- 
nificant discrepancy  from  the  data.  In  particular,  the  characteristic  flattening  of  the  measured  pressure  distri- 
bution due  to  a large  recirculation  zone  aft  of  the  comer  is  significantly  underpredicted.  Such  a “forward 
facing  step”  flow  is  well  known  to  provide  significant  challenges  to  single  phase  turbulence  models.  The  con- 
ventional model  employed  here  (high  Reynolds  number  k-e)  has  well-documented  difficulties  with  stagnated, 
high  strain  and  recirculating  flows.  All  of  these  characteristics  are  embodied  in  the  blunt  head-form  flow.  Our 
efforts  to  “remedy”  these  single  phase  turbulence  modeling  shortcomings  by  deploying  several  approaches 
that  have  appeared  in  the  literature  (q-co,  sublayer  modeling,  RNG)  have  yet  to  yield  a fully  satisfactory  reso- 
lution. As  observed  by  Shyy  [33],  [“there  are  no  quick,  practical  solutions  to  handle  this  challenge.  However, 
models  capable  of  handling  (i)  substantial  departure  from  equilibrium  between  production  and  dissipation  of 
the  turbulent  kinetic  energy,  (ii)  anisotropy  between  main  Reynolds  stress  components,  and  (iii)  turbulence- 
enhanced  mass  transfer  across  the  phase  inter-face,  should  be  emphasized”.] 

Several  parameters  of  relevance  in  the  characterization  of  cavitation  bubbles  include  body  diameter,  d,  bubble 
length,  L,  bubble  diameter,  d^,  and  form  drag  coefficient  associated  with  the  cavitator,  CD.  Some  ambiguity  is 
inherent  in  both  the  experimental  and  computational  definition  of  the  latter  three  of  these  parameters.  Bubble 
closure  location  is  difficult  to  define  due  to  unsteadiness  and  its  dependence  on  after-body  diameter  (which 
can  range  from  0 [isolated  cavitator]  to  the  cavitator  diameter).  Accordingly,  bubble  length  is  often,  and  here, 
taken  as  twice  the  distance  from  cavity  leading  edge  to  the  location  of  maximum  bubble  diameter  (see  Figure 
10).  The  form  drag  coefficient  is  taken  as  the  pressure  drag  on  an  isolated  cavitator  shape.  For  cavitators  with 
afterbodies,  such  as  here,  the  pressure  contribution  to  CD  associated  with  the  back  of  the  cavitator  is  assumed 
equal  to  the  cavity  pressure  (=  pv).  For  the  simulations,  dm  is  determined  by  examining  the  cq  = 0.5  contour 
and  determining  its  maximum  radial  location. 
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Figure  7a.  Comparison  of  predicted  and  measured  surface  Figure  7b.  Predicted  liquid  volume  fraction  and  surface 
pressure  distributions  at  several  cavitation  numbers  for  a pressure  contours,  selected  streamlines  and  computational 
hemispherical  fore-body.  grid  for  a hemispherical  fore-body,  o=0.3. 


Figure  8a.  Comparison  of  predicted  and  meas  ured  surface 
pressure  distributions  at  several  cavitation  numbers  for  a 
conical  fore-body. 


Figure  8b.  Predicted  liquid  volume  fraction  and  surface 
pressure  contours,  selected  streamlines  and  computational 
grid  for  a conical  fore-body,  0=0.3. 
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Figure  9a.  Comparison  of  predicted  and  measured  surface  Figure  9b.  Predicted  liquid  volume  fraction  and  surface 
pressure  distributions  at  several  cavitation  numbers  for  a pressure  contours,  selected  streamlines  and  computational 
blunt  fore-body.  grid  for  a blunt  fore-body,  o=0.4. 


In  Figure  11a,  the  quantity  L/(dCp  ) is  plotted  against  cavitation  number  for  a large  number  of  experimental 
data  sets  assembled  by  May  [22]  from  a variety  of  sources  and  for  numerous  simulations  made  with  six  cavi- 
tator  shapes.  That  L/(dCp  ) should  correlate  with  o has  been  long  established  theoretically  and  experimen- 
tally (Reichardt  [27],  Garabedian  [11],  for  example).  Despite  the  significant  uncertainties  associated  with 
experimental  and  computational  evaluation  of  L and  CD,  the  data  and  simulations  do  correlate  well,  close  to 
independently  of  cavitator  shape.  The  cone  shaped  cavitator  exhibits  some  underprediction  of  this  parameter 
at  higher  cavitation  numbers,  consistent  with  the  bubble  length  underprediction  observed  in  the  pressure  dis- 
tribution comparisons  presented  above. 

The  so-called  fineness  ratio  of  the  cavity,  L/ch,  is  plotted  against  cavitation  number  in  Figure  lib  for  the 
same  experimental  and  computational  data  sets.  These  parameters  again  correlate  well  for  both  experiment 
and  simulation,  though  at  higher  cavitation  numbers  there  is  a spread  in  both.  Considering  the  difficulties  in 
quantifying  smaller  bubble  sizes,  we  hesitate  to  draw  any  conclusions  as  to  the  possible  sources  of  this  spread. 


Figure  10.  Bubble  length  and  bubble  diameter  definitions. 

In  the  cases  considered  so  far,  steady-state  solutions  were  obtained  for  higher  cavitation  numbers.  At  lower 
cavitation  numbers  (i.e.  larger  bubbles),  significant  large  scale  unsteadiness  appears  in  the  aft  region  of  the 
bubble/reentrant  jet  region.  Vorticity  and  condensable  vapor  are  shed  from  this  region.  Accordingly,  transient 
simulations  (dual-time  stepping)  were  performed  when  pseudo-timestepping  failed  to  converge.  For  the  0-,  1/ 
4-,  1/2-,  1-  and  2-caliber  and  conical  cavitator  shapes,  transient  computations  were  performed  for  cavitation 
numbers  less  than  or  equal  to  o = 0.5,  0.3,  0.35,  0.15,  0.05and  0.4,  respectively.  Non-dimensional  physical 
time  steps  of  At/t^  = 0.007  were  utilized  for  these  computations.  In  Figures  11a  and  lib,  time  accurate  CFD 
results  are  represented  with  filled  symbols,  steady  state  results  with  open  symbols. 
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Ventilated  Cavitation 

Ventilated  cavities  are  also  of  interest  to  the  authors  because  for  large  enough  values  of  free  stream  pressure 
the  only  way  to  generate  large  cavities  is  to  inject  some  mixture  of  condensable  and  non-condensable  gases 
into  the  flow  near  the  body  leading  edge.  In  general,  ventilated  cavities  exhibit  similar  dynamics  to  natural 
cavities  (May  [22]).  Ventilated  cavities  do,  however,  tend  to  be  more  stable  than  natural  cavities  near  the  aft 
end  of  the  cavitation  bubble.  Also,  as  discussed  in  the  authors’  companion  paper  [38]  all  non-condensable  gas 
must  mix  with  the  free  stream  liquid  and  be  transported  downstream. 
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Figure  11a.  Comparison  of  L/(dCD  ) vs.  o for  numer- 
ous fore-body  shapes.  Experimental  data  adapted  from 
May  [22].  Open  symbols  represent  steady  computations, 
filled  symbols  represent  transient  computations. 


Figure  1 lb.  Comparison  of  L/dm  vs.  o for  numerous 
fore-body  shapes.  Experimental  data  adapted  from  May 
[22],  Open  symbols  represent  steady  computations, 
filled  symbols  represent  transient  computations. 


The  axisymmetric  blunt  fore-body  configuration  presented  above  was  run  with  no  mass  transfer  but  with  non- 
condensable gas  injection  just  aft  of  the  leading  edge.  A range  of  injection  mass  flow  rates  were  specified, 
yielding  a range  of  ventilated  cavity  sizes.  The  resulting  cavities  do  not  close  in  the  sense  that  all  vapor  is  con- 
densed, however  a distinct  bubble  shape  is  observed  whose  geometry  is  quantified  as  detailed  above.  This  is 
illustrated  in  Figure  12,  where  the  predicted  liquid  volume  fraction  fields  are  shown  for  the  blunt  head-form 
natural  and  vented  cavities  at  c = 0.4.  (The  cavitation  number  in  ventilated  cavities  is  defined  not  from  the 
vapor  pressure  but  from  the  cavity  pressure.) 

In  Figures  13a  and  1 3b,  L/(dCD  ) and  L/dm  are  plotted  for  the  blunt  head-form  at  a range  of  cavitation  num- 
bers. The  natural  cavity  results  presented  above  are  reproduced  in  these  plots  along  with  the  results  of  six  ven- 
tilated cavity  runs.  (Note  that  the  cavitation  number  is  not  explicitly  specified  in  the  ventilated  case;  rather  it  is 
an  outcome  of  specified  ventilation  flow  rate.  This  is  why  the  ventilated  data  predictions  in  Figures  13a  and 
13b  are  not  at  precisely  the  same  cavitations  numbers  as  the  natural  ventilation  cases.)  The  similarity  between 
the  natural  and  ventilated  cavity  results  are  generally  affirmed  by  the  simulations,  in  that  these  parameters  cor- 
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relate  well  with  each  other  and  with  the  experimental  data  provided  in  Figures  1 la  and  lib. 


Figure  12.  Predicted  liquid  volume  fraction  contours  for 
axisymmetric  natural  and  ventilated  cavities  about  a 
cylindrical  configuration  with  blunt  fore-body  (o  = 0.4). 
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Figure  13a.  Comparison  of  L/(dc{,  ) vs.  o for  natural  Figure  13b.  Comparison  of  L/dm  vs.  o for  natural  and 
and  ventilated  cavities  about  a blunt  fore-body.  ventilated  cavities  about  a blunt  fore-body. 


Unsteady  Characteristics  of  Cavities 

As  mentioned  above,  developed  cavities  exhibit  large  scale  unsteadiness  associated  with  re-entrant  jets,  peri- 
odic ejection  of  non-condensable  gas,  and  cavity  “pulsations”.  Accordingly,  we  have  employed  a time-accu- 
rate formulation  in  the  analysis  of  developed  cavitation.  Our  particular  interest  in  cavitator  design  has 
motivated  the  hydrodynamic  performance  parameter  model  assessments  and  validation  carried  out  above. 
There,  most  of  the  simulations  were  unsteady  and  the  presented  results  were  time  averaged. 

There  is  also  particular  interest  in  the  time  dependent  (as  opposed  to  time  averaged)  characteristics  of  devel- 
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oped  cavities,  as  these  physics  play  importantly  in  vehicle  acoustics,  body  wetting  and  bubble  cloud  collapse. 
Motivated  by  the  former  two,  we  have  carried  out  an  assessment  of  the  accuracy  of  the  predicted  temporal  sta- 
tistics provided  by  the  code. 

We  first  illustrate  the  qualitative  temporal  characteristics  of  an  analysis  of  a bluff  forebody  shape  at  moderate 
cavitation  numbers.  In  such  flows  it  is  well  known  that  the  entire  cavity  can  be  highly  unsteady,  with  “re- 
entrant” liquid  issuing  quasi-periodically  from  the  aft  end  of  the  bubble  and  traveling  all  the  way  to  the  front 
of  the  bubble.  In  Figure  14,  a time-sequence  of  predicted  vapor  volume  fraction  is  reproduced  for  a 1/4-caliber 
ogive  simulation  at  a Reynolds  number  of  1.36  x 105  and  a cavitation  number  of  0.3.  A 193  x 65  mesh  and  a 
non-dimensional  physical  time-step  of  0.007  was  used  for  this  computation.  Clearly  captured  is  the  transport 
of  a region  of  liquid  towards  the  front  of  the  cavity.  There,  the  liquid  interacts  with  the  bubble  leading  edge, 
the  top  of  this  liquid  region  being  sheared  aftward  while  the  bulk  of  the  fluid  proceeds  upstream  “pinching 
off’  the  bubble  near  the  leading  edge. 


Figure  14.  Time  sequence  of  predicted  vapor  volume 
fraction  for  flow  over  a 1/4  caliber  ogive  with  cylindri- 
cal afterbody,  a = 0.3 
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Figure  15.  a)  Pseudo-time  convergence  history  at  four 
successive  physical  time  steps  for  transient  flow  over  a 
1/4  caliber  ogive  with  cylindrical  afterbody,  a = 0.3. 
b)  Segment  of  time  history  of  predicted  CD  for  this  case. 


Figure  15  shows  two  other  elements  of  this  particular  transient  simulation.  The  inner-  or  pseudo-time  conver- 
gence history  for  four  successive  physical  time  steps  is  shown  in  Figure  15a.  Nearly  a three  order-of-magni- 
tude  drop  in  the  axial  velocity  residuals  is  obtained  at  each  physical  time  step  using  10  pseudo-time-steps. 
Figure  1 5b  shows  a segment  of  the  time  history  of  predicted  drag  coefficient  for  this  case 

Stinebring  et  al.  [36]  documented  the  unsteady  cycling  behavior  of  several  axisymmetric  cavitators.  Their 
report  included  results  for  both  ventilated  and  natural  cavitation.  The  unsteady  performance  of  45°  (22.5° 
half-angle)  conical,  hemispherical,  and  0-caliber  ogival  cavitators  at  a range  of  cavitation  numbers  were  docu- 
mented. Natural  cavitation  analysis  comparisons  have  been  included  here. 

Figure  16  contains  a series  of  snapshots  of  the  predicted  volume  fraction  field  from  an  unsteady  model  com- 
putation of  flow  over  a 0-caliber  (blunt)  cavitator.  Flere  the  Reynolds  number  (based  on  diameter)  was 
1.46xl05  and  the  cavitation  number  was  0.3.  This  result  is  presented  over  an  approximate  model  cycle.  The 
figure  also  includes  the  corresponding  time  segment  of  drag  coefficient.  Note  that  the  spikes  in  drag  near 
t=37.725  and  t=38.925  seconds  correspond  to  reductions  in  the  relative  amount  of  vapor  near  the  sharp  lead- 
ing edge.  This  marks  the  progress  of  a bulk  volume  of  liquid  from  the  closure  region  to  the  forward  end  of  the 
cavity  as  part  of  the  reentrant  jet  process.  Although  far  from  regular,  these  spikes  also  delineate  the  approxi- 
mate model  cycle. 

The  computed  physics  in  Figures  14  and  23  correspond  qualitatively  to  film  footage  of  blunt  cavitators  at 
intermediate  cavitation  numbers.  This  is  illustrated  in  Figure  17  which  is  a photograph  (adapted  from  Stine- 
bring [37])  of  a 0-caliber  axisymmetric  cavitator  operating  at  ReD=2.9xl05,  Gs0.35.  This  picture  serves  to 
illustrate  the  basic  phenomenon  of  natural  sheet  cavitation  as  it  is  captured  by  the  model.  This  result  is  notable 
for  the  spatially  and  temporally  irregular  nature  of  the  computed  flow  field.  Even  after  significant  integration 
effort,  a clearly  periodic  result  had  not  emerged.  Thus,  to  deduce  the  dominant  frequency  with  some  confi- 
dence, it  was  necessary  to  apply  ensemble  averaging. 

Note,  in  Figure  16  that  over  a significant  portion  of  the  sequence,  the  leading,  or  formative,  edge  of  the  cavity 
sits  slightly  downstream  from  and  not  attached  to  the  sharp  corner.  In  their  experiments,  Rouse  and  McNown 
observed  this  phenomenon.  They  suggested  that  this  delay  in  cavity  formation  was  due  to  the  tight  separation 
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eddy  which  forms  immediately  downstream  of  the  corner  and,  hence,  locally  increases  the  pressure.  The  cor- 
responding evolution  of  cavitation  further  downstream,  at  the  separation  interface,  was  proposed  to  be  due  to 
tiny  vortices.  These  vortices,  after  some  time,  subsequently  initiate  the  cavity. 

Figure  18  shows  a single  frame  at  t=37.8  seconds  from  the  same  model  calculation  (as  shown  in  Figure  23). 
Here,  to  clarify  what  is  captured,  the  volume  fraction  contours  have  been  enhanced  with  illustrative  stream- 
lines. Note  that  these  are  streamlines  drawn  from  a frozen  time  slice.  Nonetheless,  if  all  of  the  details  envi- 
sioned by  Rouse  and  McNown  were  present,  the  streamlines  should  indicate  smaller/tighter  vortical  flows. 
The  current  level  of  modeling  was  unable  to  capture  small  vortical  structures  in  the  flow.  However,  the  overall 
computation  was  apparently  able  to  capture  the  gross  affects  of  these  phenomena  and  reproduce  a delayed 
cavity.  In  fact  from  examination  of  the  cavity  cycle  evolution  shown  in  Figure  23,  and  the  streamlines  shown 
in  the  snapshot,  it  appears  that  gross  unsteadiness  is  driven  by  a combination  of  a reentrant  jet  and  some  type 
of  cavity  pinching.  The  pinching  process  is  particularly  well  demonstrated  in  Figure  23  from  t=3 8.125  to 
38.325  seconds.  However,  rather  than  complete  division  and  convection  into  the  free  stream,  it  should  be 
noted  that,  in  later  frames  of  Figure  23,  the  pinched  portion  of  the  cavity  appears  to  rejoin  the  main  cavity 
region. 

The  low  frequency  mode  apparent  in  most  of  the  experimental  0-caliber  results  appears  to  have  been  captured 
at  the  lowest  cavitation  number  la  0.3),  as  shown  in  Figure  23,  and  is  evidenced  in  the  test  photograph  (Fig- 
ure 17).  In  Figure  19,  the  drag  coefficient  history  for  a 40  model  second  interval  from  the  same  computation  as 
in  Figure  23  is  shown.  Here,  a clear  picture  of  the  persistence,  over  a long  integration  time,  of  the  irregular 
flow  behavior  is  documented.  At  higher  cavitation  numbers,  the  current  set  of  0-caliber  cavitator  results  indi- 
cate a more  regular  periodic  motion.  This  is  contrary  to  the  experimental  data.  However,  as  Figure  18  indi- 
cates, the  ability  to  capture  this  motion  at  any  cavitation  number  may  not  necessarily  require  the  explicit 
capture  of  the  finer  flow  details  of  the  vortical  flow  structure.  This  is  encouraging  and  suggests  that  with 
increased  computational  effort,  without  altering  the  current  physical  model,  the  representation  of  this  phenom- 
enon could  be  improved  over  a greater  range  of  cavitation  numbers. 

Figure  20  presents  the  spectral  content  of  the  result  given  in  Figure  19.  This  power  spectral  density  plot  is 
based  on  four  averaged  Hanning  windowed  data  blocks  of  the  time  domain  result.  To  eliminate  the  start-up 
transient  effect,  the  record  was  truncated,  starting  at  t=10  seconds  and,  to  tighten  the  resulting  confidence 
intervals,  more  time  domain  results,  after  t=40  seconds  were  included.  As  is  typical  of  highly  nonlinear 
sequences,  the  experience  of  this  unsteady  time  integration  demonstrated  that,  additional  time  records  merely 
enrich  the  power  spectral  density  function.  However,  the  additional  records  do  serve  to  improve  the  confi- 
dence intervals,  and,  therefore,  add  reliability  to  the  numerical  convergence  process.  The  model  result  used 
was,  as  indicated  by  the  confidence  intervals,  sufficient  for  a comparison  to  experimental,  unsteady  results. 
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Figure  16.  Modeled  flow  over  a 0-caliber  ogive.  Liq- 
uid volume  fraction  contours  and  corresponding  drag 
history.  o=03.  ReD=1.46xl05. 


Figure  17.  Zero  caliber  ogive  in  water  tunnel  at 
ReD=2.9xl05,  os0.35  (adapted  from  Stinebring  [37]). 
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Figure  18.  Snapshot  of  modeled  flow  over  a 0-caliber 
ogive.  Liquid  volume  fraction  contours  and  selected 
streamlines.  c=0.3.  Reo=1.46xl05. 


Figure  19.  Model  time  record  of  drag  coefficient  for  flow  over  a 0-caliber  ogive  at  ReD=l  .46xl05  and  a=0.3.  In  model 

units,  D/Um  = 0.146  (s),  physical  time  step,  At  = 0.001  (s). 


Figure  20.  Results  for  0-caliber  ogive  at  ReD=1.46xl05  and  o= 0.3.  Power  spectral  density  function  with  50%  confidence 

intervals  shown. 


Figure  21  contains  a series  of  snapshots  from  the  unsteady  model  computation  of  a hemispherical  cavitator  at 
a Reynolds  number  (based  on  diameter)  of  1.36x10''  and  a cavitation  number  of  0.2.  This  result  is  presented 
over  a period  slightly  longer  than  the  approximate  model  cycle.  In  this  case  the  model  Strouhal  frequency  is 
0.0326.  There  are  ten  frames  presented,  and  the  first  (or  last)  nine  of  those  ten  constitute  an  approximate 
model  cycle.  The  drag  history  trace  in  Figure  22  demonstrates  how,  relative  to  the  modeled  flow  over  the 
blunt  forebody,  the  pattern  of  flow  over  the  hemispherical  forebody  is  regular  and  periodic.  This  is  consistent 
with  experimental  observations  made  (for  example)  by  Rouse  and  Mcnown  [29],  Note  the  evolution  of  flow 
shown  in  Figure  21  as  it  compares  to  the  drag  history  shown  in  Figure  22.  As  would  be  expected,  the  large 
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spike  in  drag  corresponds  to  the  minimum  in  vapor  shown  near  the  modeled  t=1.6  seconds. 

Figure  23  contains  a time  record  of  drag  coefficient  during  modeled  flow  over  a conical  forebody  and  cylinder 
at  a Reynolds  number  of  1.36x10'''  and  cavitation  number  of  0.2.  The  Strouhal  frequency  based  on  this  result  is 
0.0383.  As  anticipated,  due  to  the  expected  stability  of  cavities  about  this  shape,  this  model  flow  exhibited 
very  regular  cycling  with  little  additional  strong  components  from  secondary  modes. 


Figure  2 1 . Liquid  volume  fraction  countours.  Modeled  flow  over  a hemispherical  forebody 

and  cylinder.  o=0.2,  ReD=1.36xl05. 


Figure  24  contains  a large  survey  of  unsteady  computational  and  experimentally  obtained  data  [36].  The 
numerical  results  in  this  figure  summarize  this  validation  effort.  Here,  Strouhal  frequency  is  shown  over  a 
range  of  cavitation  numbers.  Computational  results  are  given  for  hemispherical,  1/4-caliber,  conical,  and  0- 
caliber  forebodies.  Unsteady  experimental  data  is  included  for  the  hemispherical,  conical  and  0-caliber 
shapes.  Computational  results  for  the  hemisphere,  1/4-caliber  and  conical  forebodies,  were  obtained  at  a Rey- 
nolds number  based  on  diameter  of  1.36xl05.  For  the  0-caliber  ogive,  computations  were  made  at  a Reynolds 
number  of  1.46xl05.  In  addition,  for  the  hemisphere,  results  are  included  for  Reynolds  numbers  of  1.36xl06 

■7 

and  1.36x10  . The  experimental  results  included  in  the  figure  were  obtained  at  Reynolds  numbers  ranging 
from  3.5x10s  to  1.55xl06. 
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Figure  22.  Unsteady  drag  coefficient.  Flow  over  a 
hemispherical  forebody  and  cylinder.  o=0.2, 
ReD=1.36xl05.  In  model  units,  D/U^  = 0.136(s), 
physical  time  step.  At  = 0.001  (s). 


Figure  23.  Predicted  time  record  of  drag  coefficient 
for  flow  over  a conical  forebody  and  cylinder  at 
Re[)=1.36xl05  and  o=0.3.  In  model  units, 

D/Um  = 0.136  (s),  physical  time  step, 

At  = 0.0025 

As  a blanket  observation,  the  spread  of  data  between  the  experiments  and  computations  in  Figure  24  is  signif- 
icant. However,  there  are  several  encouraging  items  to  be  reviewed.  It  is  clear  that  (for  a given  cavitation  num- 
ber) the  computational  results  are  bounded  by  the  experimental  data,  and  the  proper  trends  (rate  of  change  of 
Strouhal  frequency  with  cavitation  number)  are  well  captured.  More  insight  into  the  physical  relevance  of  the 
data  requires  examination  of  specific  results. 

For  the  hemispherical  forebody  results,  as  may  be  seen  in  Figure  24,  there  is  a significant  but  almost  constant 
offset  between  the  measured  unsteady  data  and  the  modeled  results  both  of  which  appear  to  follow  a linear 
trend  over  the  range  presented.  An  interesting  result  occurs  in  the  model  data  for  the  hemispherical  forebody 
with  a Reynolds  number  of  1.36x10  (pentagrams  in  Figure  24).  Here  the  numerical  results  appear  to  agree 
quite  well  with  the  experimental  data  for  hemispherical  forebodies.  The  experiments  were  taken  at  an  order  of 
magnitude  lower  Reynolds  number,  but  the  agreement  is  apparent  in  both  cases  where  model  results  have 
been  obtained.  For  design  purposes,  this  may  suggest  an  avenue  towards  model  calibration. 
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Figure  24.  Axisymmetric  running  cavitators  with  cylindrical  afterbodies.  Strouhal  frequency  and  cavitation  number. 
Model  results  (open  symbols)  and  data  reported  in  Stinebring  36. 


Another  result  found  in  the  Str  versus  0 plot  (Figure  24)  is  the  tendency  of  the  modeled  flows  to  become 
steady  at  higher  cavitation  numbers.  For  the  0-caliber  or  the  conical  cavitators,  this  is  the  reason  model  results 
are  not  included  for  cavitation  numbers  greater  than  0.4.  For  the  modeled  hemisphere,  the  upper  limit  of  cavi- 
tation number  to  yield  unsteady  model  results  was  found  to  be  Reynolds  number  dependent.  At  a 
ReD=1.36xl05,  the  maximum  cavitation  number  yielding  an  unsteady  result  was  a = 0.35,  at  ReD=1.36xl06, 
that  number  was  o = 0.45,  and  at  ReD=1.36xl07,  the  maximum  cavitation  number  for  unsteady  computations 
was  not  determined.  This  result  may  indicate  a limit  of  the  computational  grid  applied  to  the  problems  rather 
than  a limit  of  the  level  of  physical  modeling.  In  addition,  physically  in  the  mode  of  unsteadiness  present,  a 
transition  does  occur  from  cavity  driven  to  separated,  turbulent,  but  single  phase  driven  flow. 

For  the  conical  forebody,  the  datum  shown  in  Figure  24  suggests  that  the  cycling  frequency  should  be  higher, 
0.123.  It  is  worth  considering  that  the  Reynolds  number  of  the  experimental  flow  was  3.9xl05  and  that  the 
general  trend  with  increasing  Reynolds  numbers  is  to  increased  frequency.  However,  based  on  the  standard 
level  of  dependence  of  Strouhal  frequency  (see  Schlichting  [31]  for  example)  on  Reynolds  number  for  bluff 
body  flows,  it  would  seem  unlikely  that  the  rate  of  change  in  frequency  with  Reynolds  number  (at  ReD  = 105 ) 
would  be  as  high  as  three  to  two.  In  addition,  compared  to  shapes  with  geometrically  smooth  surfaces,  the 
nature  of  unsteady  flow  over  a conical  shape  is  not  expected  to  be  nearly  so  dependent  on  Reynolds  number. 
In  the  case  of  a cone,  at  low  values  of  cavitation  number  (i.e.  0=0.3),  the  separation  location,  and,  hence,  the 
likely  forward  location  of  the  cavity,  is  rarely  in  question. 

A trend  that  is  captured  in  the  model  results  but  not  represented  in  the  experimental  data  included  here,  is  the 
tendency  for  the  Strouhal  frequency  of  a given  cavitator  shape  to  exhibit  two  distinct  flow  regimes.  The  first 
regime  exists  at  moderate  cavitation  numbers  and  is  indicated  by  a low  Strouhal  frequency  where  the  value  of 
Str  will  have  an  apparent  linear  dependence  on  0.  The  second  regime  tends  toward  much  higher  cycling  fre- 
quencies. Here  the  dependent  Strouhal  frequency  appears  to  asymptotically  approach  a vertical  line  with 
higher  cavitation  number,  just  prior  to  the  complete  elimination  of  the  cavity.  This  is  documented  in  Stine- 
bring [36]  and  demonstrated  in  Figure  24  for  the  modeled  hemisphere  at  ReD=1.36xl06.  Based  on  the  model 
results,  it  appears  that  this  is  characteristic  of  a change  from  a flow  mode  dominated  by  a large  unsteady  cav- 
ity to  one  dominated  by  other,  single-phase,  turbulent,  sources  of  unsteadiness. 

During  this  investigation,  some  effort  towards  the  establishment  of  temporal  and  spatial  discretization  inde- 
pendence was  made.  As  a requirement  of  the  model,  to  accommodate  the  use  of  wall  functions,  for  regions  of 
attached  liquid  flow,  fine-grid  near-wall  points  were  established  at  locations  yielding  10<y+<100. 

Temporal  convergence  was  established  by  the  successive  reduction  of  time  integration  step  for  several 
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selected  cases.  Figure  25  contains  a comparison  of  the  spectral  content  of  results  for  flow  over  a hemispherical 
forebody  and  cylindrical  afterbody,  with  ReD=1.36xl05  and  0=0.3,  for  three,  successively  smaller,  integration 
step  sizes.  Here,  with  a physical  time  step,  At  = 0.005  seconds,  the  computation  resulted  in  a Strouhal  fre- 
quency, Str=0.0680,  with  a time  step,  At  = 0.0025  seconds,  Str=0.0622,  and  with  At  = 0.001  seconds, 
Str=0.0680.  More  significantly,  as  demonstrated  in  the  figure,  for  the  smaller  two  integration  step  sizes,  over 
the  range  of  relevant  (shown)  frequencies,  there  was  very  similar  modal  behavior. 

Unfortunately  only  the  fine-grid  models  tended  to  provide  unsteady  results.  Thus  time  and  spatial  fidelity 
were  judged  independently.  A demonstration  of  the  steady-state  spatial  convergence  of  the  modeled  conical 
forebody  and  cylindrical  afterbody  is  given  in  Figure  5. 


Figure  25.  Spectral  comparison  of  effect  of  physical  integration  time  step  size  on  Q history.  Flow  over  a hemispherical 

forebody  with  cylindrical  afterbody.  ReD=T.36xl05.  o= 0.3. 


Three-Dimensional  Analyses 

Natural  Cavitation  about  1/2-  and  1-Caliber  Ogives 

In  order  to  demonstrate  the  three-dimensional  capability  of  the  method,  a model  of  the  hemispherical  fore- 
body configuration  studied  above  was  run  at  numerous  angles  of  attack  and  a cavitation  number  of  0.3.  A 
97x33x65  mesh  was  utilized  (corresponding  to  the  “medium”  mesh  size  discussed  in  grid  studies  above).  The 
domain  was  decomposed  into  8 subdomains  azimuthally  and  run  on  8 SGI  RSI  OK  Octane  machines.  Parallel 
efficiencies  of  85%  were  achieved  for  these  problems. 

Figure  26  provides  sample  results  for  angles  of  attack  of  0.0°,  2.5°,  5.0°  and  7.5°.  These  plots  include  pressure 
contours  on  the  plane  of  symmetry,  sample  streamlines  and  the  cavitation  bubble  shape  as  identified  with  an 
isosurface  of  a/  = 0.99.  Several  interesting  features  are  observed  in  the  predictions.  In  particular,  the  flows  are 
seen  to  be  highly  three-dimensional  in  nature  at  angle-of-attack.  A recirculation  zone  aft  of  the  bubble,  grows 
with  angle  of  attack.  This  diminishes  the  local  pressure  recovery  associated  with  the  bubble-induced  blockage 
and  this  in  turn  leads  to  a local  collapse  of  the  bubble  on  the  top  of  the  body.  Indeed  at  angle-of-attack  the  bub- 
ble is  seen  to  have  its  greatest  axial  extent  off  of  the  symmetry  plane  of  the  geometry. 

Figure  27  shows  the  results  of  a three-dimensional  simulation  of  cavitating  flow  over  a 1-caliber  ogive  fore- 
body with  cylindrical  afterbody  at  a 10°  angle-of-attack.  A cavitation  number  of  0.32  was  specified.  Again  a 
97  x 33  x 65  mesh  was  utilized,  the  domain  was  decomposed  into  eight  subdomains  azimuthally  and  run  on 
eight  processors  and  at  this  cavitation  number  using  this  grid  a steady  state  solution  could  be  obtained.  Figure 
27  illustrates  the  predicted  bubble  shape  and  streamline  pattern  for  this  simulation.  The  bubble  shape  is  highly 
three-dimensional  in  nature;  it  does  not  close  on  the  pressure  side  of  the  body.  The  streamline  pattern  is  char- 
acterized by  a large  recirculation  zone  aft  of  the  bubble  and  significant  azimuthally  oriented  vortical  struc- 
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tures. 

Figure  28  shows  the  convergence  history  for  this  case.  A three  order-of-magnitude  drop  in  the  axial  velocity 
residual  is  achieved  in  800  pseudo-time-steps.  It  is  observed  that  the  8-block/8-processor  simulation  exhibits 
nearly  identical  convergence  behavior  to  a 1 -block/ 1 -processor  run. 

Interaction  of  a Control  Surface  with  Phase-Separated  Non-condensable  Gas  and  Liquid  Impinge- 
ment Streams 

The  second  three-dimensional  simulation  presented  is  that  of  a wedge  shaped  control  surface  interacting  with 
an  incoming  stream  of  phase-separated  water  and  air.  This  configuration  was  tested  by  the  fourth  author  in  the 
12”  water  tunnel  at  the  Penn  State  Applied  Research  Laboratory  (unpublished).  The  test  was  run  with  co- 
directed air  and  water  streams.  Air  was  injected  along  the  bottom  of  the  tunnel  such  that  the  unperturbed  gas- 
liquid  interface  impinged  upon  the  sharp  leading  edge  at  about  25%  span  from  the  base  of  the  fin.  Air  and 
water  velocities  were  approximately  the  same.  The  configuration  was  tested  at  a range  of  angles-of-attack. 

A 189,546  vertex  grid  was  used  for  the  results  presented  here.  The  simulation  was  run  on  16  processors  with  a 
10°  angle-of-attack  and  a cavitation  number  of  0.15.  Using  this  relatively  coarse  mesh  and  the  k-e  model,  a 
steady-state  solution  was  obtained,  despite  the  presence  of  a blunt  trailing  edge  and  its  associated  recirculation 
zone  along  the  span.  As  we  refine  our  analyses  of  this  class  of  application,  time  accurate  simulations  using 
upwards  of  106  nodes  will  be  required. 

These  flows  are  characterized  by  several  physical  features  of  academic  and  practical  interest.  At  angle-of- 
attack,  the  liquid  and  gas  streams  are  turned  through  approximately  the  same  angle,  but  because  the  density 
ratio  is  high  (=1000  here),  significant  spanwise  pressure  gradients  arise  along  the  lifting  surfaces.  Accord- 
ingly, on  the  pressure  side,  the  gas-liquid  interface  is  deflected  downward,  giving  rise  to  deceleration  and 
acceleration  of  the  liquid  and  gas  streams  respectively.  On  the  suction  surface  the  lower  pressure  in  the  liquid 
gives  rise  to  an  upward  deflection  of  the  gas  interface.  Attendant  to  these  interface  deflections  is  a loss  in  lift. 

Adjacent  to  the  sharp  leading  edge  on  the  suction  side,  the  local  static  pressure  becomes  low  due  to  leading 
edge  separation.  Local  natural  cavitation  occurs,  and  this  vapor  merges  with  the  swept  up  non-condensable 
interface  so  that  the  entire  suction  side  is  enveloped  in  gas  phase.  Aft  of  the  blunt  trailing  edge,  the  flow  is 
recirculating  and  therefore  the  local  static  pressure  is  also  low.  This  gives  rise  to  a sweeping  up  of  the  pressure 
side  cavity  and  some  natural  cavitation  so  that  the  entire  wake  region  is  principally  gas  phase. 

Figure  29  shows  an  oblique  front/top  view  of  the  predicted  flow  field.  Each  of  the  flow  features  discussed 
above  is  clearly  observed  in  the  simulation.  Of  particular  note  is  that  the  three-species  formulation  enables  the 
separate  prediction  and  identification  of  vaporous  and  non-condensable  regions  of  the  flow.  Figure  30a  shows 
a suction  side  video  frame  of  the  tested  flow  configuration  at  this  angle-of-attack.  The  enveloping  of  the  suc- 
tion surface  in  gas  phase  and  the  primarily  gas  phase  wake  are  observed.  Figure  30b  shows  a pressure  side 
photograph.  A similar  view  of  the  simulation  is  presented  in  Figure  30c  which  shows  that  the  pressure  surface 
interface  deflection  and  gas-vapor  wake  region  are  qualitatively  consistent  with  experimental  observation. 


Figure  26.  Predicted  3-D  flow  field  with  natural  cavitation  about  hemispherical  fore-body  at  several  angles  of 
attack.  0 = 0.3.  Liquid  volume  fraction  = 0.99  isosurface. 


single  block 
multiblock 


Iteration 


Figure  27.  Predicted  bubble  shape  (designated  by  a;  = 0.99 
isosurface)  and  streamlines  for  a 1 -caliber  ogive  at  a 1 0° 
angle-of-attack. 


Figure  28.  Comparison  of  1 -block/ 1 -processor 
and  8-block/8-processor  convergence  histories  for 
a 1 -caliber  ogive  at  a 10°  angle-of-attack. 
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Figure  29.  Front/top  view  of  CFD  simulation  of  super- 
cavitating  fin  configuration  at  1 0°  angle-of-attack. 

Isocontours  of  ang  = 0.9,  av  = 0.9  designate  non-con- 
densable gas  and  vapor  regions  respectively. 

Three-dimensional  Flow  over  a Vehicle  with  Ventilation  and  Natural  Cavitation 

Practical  high  speed  supercavitating  vehicle  designs  employ  ventilation  at  or  near  the  cavitator  leading  edge, 
where  natural  cavitation  will  likely  also  occur.  Figure  31  illustrates  elements  of  a 48-processor  three-dimen- 
sional analysis  of  a notional  vehicle  with  natural  cavitation,  noncondensable  gas  injection  and  propulsion 
stream.  This  illustrates  the  three-species  nature  of  the  formulation. 

Unsteady  Three-Dimensional  Flow  Over  a 0-Caliber  Ogive 

As  pointed  out  by  Edwards  [10],  a potential  cause  of  the  discrepancies  observed  between  predicted  and  mea- 
sured Strouhal  numbers  for  the  cavitator  series  simulations  is  the  axisymmetric  nature  of  the  unsteady  runs. 
We  have  begun  to  investigate  this  by  simulating  these  flows  three-dimensionally.  Figure  32  illustrates  ele- 
ments of  a 1,245,184  cell,  78  processor  unsteady  simulation  of  flow  over  a blunt  ogive  at  o=0.30. 
ReD=1.46xl05.  Clearly  observable  are  non-axisymmetric  predicted  cavity  shapes.  We  continue  to  investigate 
these  modes  in  the  hope  of  quantifying  them. 
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Figure  30.  a)  Suction  side  video  frame  of  super-cavitating  fin  at 
10°  angle-of-attack.  b)  Pressure  side  photograph,  c)  Pressure  side 
view  of  CFD  simulation. 
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Figure  31.  Elements  of  three-dimensional  supercavitat- 
ing  vehicle  simulations  with  ventilated  and  natural  cavi- 
tation. a)  liquid  volume  fraction  contours  and  stream 
lines  (a  = 0°),  b)  view  of  3D  geometry  and  grid,  c)  liq- 
uid volume  fraction  = 0.9  isovolume  and  vertical  veloc- 
ity contour  plane  (a  = 5°). 


Figure  32.  Three-dimensional,  turbulent,  unsteady,  two-phase  result.  1 ,245, 1 84  cell  grid.  Flow  over  blunt  ogive  with 
an  isosurface  of  volume  fraction,  a/=0.9,  and  selected  streamlines  at  two  timesteps.  <5=0.30.  ReD=1.46xl05. 

Sheet  Cavitation  in  a Centrifugal  Pump 

The  authors  have  recently  extended  our  efforts  towards  the  analysis  of  sheet  cavitation  in  turbomachinery. 
Some  representative  results  are  presented  here.  Figure  33  shows  three  simulations  of  a generic  backswept 
impeller  at  three  cavitation  numbers.  These  plots  show  a front  view  of  the  impeller,  with  surface  pressure  con- 
tours on  the  blade  surfaces  and  hub.  Also  shown  are  the  predicted  sheet  cavities  on  the  blade  surfaces  as  des- 
ignated with  an  isovolume  contour  of  liquid  volume  fraction  = 0.9.  A single  blade  passage  was  computed;  the 
results  have  been  copied  several  times  to  yield  a complete  annulus  view. 

We  have  not  yet  analyzed  the  cavitation  performance  of  this  three  dimensional  configuration.  This  is  because 
when  employing  the  higher  order  convection  differencing  and  fine  grids  required  for  accurate  cavitation  bub- 
ble size  predictions  for  a specified  cavitation  number,  the  cavities  and  thereby  the  simulation  are  unsteady.  As 
a step  towards  the  goal  of  predicting  cavitation  performance  in  pumps,  we  have  analyzed  a quasi-three-dimen- 
sional midspan  streamsheet  representation  of  the  same  configuration.  Performance  predictions  for  this  series 
of  analyses  are  shown  in  Figure  34.  In  Figure  34a,  head  coefficient  versus  flow  coefficient  is  plotted  for  a 
spectrum  of  flow  coefficients  and  cavitation  numbers.  Single  phase  results  are  plotted  as  solid  red  diamonds. 
A characteristic  head  drop  off  is  observed  off  design  (though  clearly  efficiency  levels  and  off  design  trends  are 
unrealistic  due  to  the  quasi-three-dimensional  modeling  employed).  A sequence  of  successively  lower  cavita- 
tion numbers  were  run  for  several  flow  coefficients.  The  reduction  in  predicted  head  attendant  to  the  sheet 
cavitation  are  also  plotted  in  Figure  34a.  Figure  34b  illustrates  the  ability  of  the  analysis  to  capture  incipient 
cavitation  breakdown  physics.  There,  head  coefficient  is  plotted  versus  cavitation  number  for  the  series  of 
quasi-three-dimensional  simulations.  At  each  flow  coefficient,  a rapid  drop-off  in  head  is  observed  at  a critical 
cavitation  number. 

The  authors  point  out  that  these  pump  cavitation  results  are  preliminary  and  several  issues  remain  to  be 
resolved  including  the  need  for  improved  mass  transfer  modeling  to  account  for  thermal  effects  on  cavitation 
breakdown,  and  three-dimensional  performance  analysis. 
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Figure  33.  Results  of  computation  of  a generic  back- 
swept impeller  at  three  cavitation  numbers.  Surface  pres- 
sure contours  on  the  blade  surfaces  and  hub  and 
isovolume  contour  of  liquid  volume  fraction  = 0.9. 


a=0.29 
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Figure  34.  Results  of  quasi-three-dimensional  midspan  streamsheet  analyses  of  a generic  backswept  impeller,  a)  pre- 
dicted efficiency  vs.  flow  coefficient,  b)  predicted  head  coefficient  vs.  cavitation  number  at  several  flow  coefficients. 


Compressible  Simulations 

Some  recently  obtained  compressible  flow  results  are  presented  in  Figures  35  and  36.  Figure  35  contains  an 
illustrative  result  of  a projectile  traveling  supersonically  relative  to  the  speed  of  sound  in  water.  On  the  left  is 
a photograph  of  the  actual  test.  At  right  is  the  model  result  for  this  case.  In  the  model  result,  the  liquid  to  vapor 
density  ratio  is  nominally  1000,  and  the  Mach  Number  is  1.03.  The  example  is  illustrative  of  challenging 
modeling  issues.  In  the  liquid  flow,  there  are  shock  waves.  These  are  resolved  with  MUSCL  interpolation  and 
the  van  Albada  limiter.  In  addition,  due  to  the  high  flow  velocity,  the  cavitation  index  is  approximately  10~4. 
In  the  figure,  dark  blue  indicates  density  of  less  than  one  and  red  indicates  density  of  nominally  1000.  Thus, 
most  of  the  flow  immediately  adjacent  to  the  object  is  completely  vaporized  as  is  the  downstream  wake.  This 
is  exhibited  in  both  the  photograph  and  the  model  result. 

Figure  36  contains  a modeled  axisymmetric  plume  from  a hypothetical  rocket  operating  underwater.  Here  the 
jet  is  supersonic  and  slightly  underexpanded.  It  is  surrounded  by  a second  gas  stream  which  is  subsonic 
(notionally  the  vehicle  cavity)  and  finally  by  a subsonic  water  free  stream.  In  this  case  the  nominal  liquid  to 
gas  density  ratio  is  again  1000.  In  the  density  field,  red  indicates  a density  of  nominally  1000,  and  blue  indi- 
cates a density  of  one  or  less.  In  the  shock  function  field,  the  classic  expansion  pattern  followed  by  shocks  are 
exhibited.  The  interaction  of  the  compressible  gas  stream  with  the  nearly  incompressible  liquid  is  demon- 
strated by  the  contraction  and  expansion,  from  top  to  bottom  in  the  figure,  of  the  gas  stream.  In  addition,  the 
interface  between  the  gas  and  liquid  contains  fully  supersonic  flow.  This  region  is  highlighted  by  the  shock 
function  due  to  the  density  gradient,  but  does  not  necessarily  contain  shocks.  A line  segment  drawn  locally 
normal  to  the  interface  would  be  nearly  isobaric.  Further  results  of  our  compressible  multiphase  efforts  are 
available  in  [20]  and  [41], 
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Figure  35.  Hypervelocity  projectile  test  and  computation,  a)  Photograph  of  a projectile  at  M=1 .03  with  respect  to  water 
sound  speed,  b)  Corresponding  simulation  showing  surface  pressure  and  field  mixture  density  contours. 

Conclusions 


A multi-phase  CFD  method  has  been  presented  and  applied  to  a number  of  high  density  ratio  developed-  and 
super-cavitating  flows.  Several  aspects  of  the  method  were  outlined  and  demonstrated  that  enable  convergent, 
accurate  and  efficient  simulations  of  these  flows.  These  include  a differential  model  and  preconditioning  strat- 
egy with  favorable  eigensystem  characteristics,  a block  implicit  dual-time  solution  strategy,  a three  species 
formulation  that  separately  accounts  for  condensable  and  non-condensable  gases,  higher  order  flux  differenc- 
ing with  limiters  and  the  embedding  of  this  scheme  in  a parallel  multi-block  Navier-Stokes  platform. 

The  two-dimensional/axisymmetric  simulations  presented  verify  the  ability  of  the  tool  to  accurately  analyze 
steady-state  and  transient  sheet-  and  super-cavity  flows.  The  three-dimensional  and  compressible  capability  of 
the  code  was  demonstrated  as  well. 

As  the  authors  proceed  with  this  research,  we  are  focusing  on  several  areas  including:  1)  validation  of  the 
method  for  compressible  constituent  fields,  2)  improved  physical  models  for  mass  transfer  and  turbulence,  3) 
extended  application  and  validation  for  steady  and  transient  three-dimensional  flows  and  4)  improved  error 
damping  through  preconditioning  and  pseudo-time-stepping  formulations  that  locally  adapt  to  problem 
parameters. 
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Figure  36.  Cartoon  vehicle  and  3-stream,  axisymmetric  plume  computation.  Center  jet  (diameter=l,  M=3  with 
respect  to  propellant  sound  speed)  surrounded  by  cavity  gas  at  free  stream  velocity  (outer  diameter=2,  M«l),  sur- 
rounded by  freestream  liquid  (Liquid  to  gas  density  ratio  1 000,  M«  1 ).  a)  mixture  density  contours,  b)  shock  func- 
tion, V Vp  , contours. 
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